Distributed self-regulation of living tissue. Effects of nonideality 
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Self-regulation of living tissue as an example of self-organization phenomena in active fractal 
systems of biological, ecological, and social nature is under consideration. The characteristic feature 
of these systems is the absence of any governing center and, thereby, their self-regulation is based 
on a cooperative interaction of all the elements. The paper develops a mathematical theory of a 
vascular network response to local effects on scales of individual units of peripheral circulation. 

First, it formulates a model for the self-processing of information about the cellular tissue state 
and cooperative interaction of blood vessels governing redistribution of blood flow over the vascular 
network. Mass conservation (conservation of blood flow as well as transported biochemical com- 
pounds) plays the key role in implementing these processes. The vascular network is considered to 
be of the tree form and the blood vessels are assumed to respond individually to an activator in 
blood flowing though them. 

Second, the constructed governing equations are analyzed numerically. It is shown that at the 
first approximation the blood perfusion rate depends locally on the activator concentration in the 
cellular tissue, which is due to the hierarchical structure of the vascular network. Then the distinction 
between the reaction threshold of individual vessels and that of the vascular network as a whole is 
demonstrated. In addition, the nonlocal component of the dependence of the blood perfusion rate 
on the activator concentration is found to change its form as the activator concentration increases. 



I. DISTRIBUTED SELF-REGULATION IN 
ACTIVE HIERARCHICAL SYSTEMS 

For a wide class of biological, ecological and social sys- 
tems, including economic ones a generalized notion of 
"homeostasis" can be introduced. By this term we under- 
stand the internal conditions required for such a system 
to survive as well as its ability to maintain them. Living 
multicellular organisms can exist if only the tempera- 
ture, oxygen concentration, etc. are in certain intervals, 
giving rise to a large number of mechanisms controlling 
these conditions (for introduction see, e.g., [l.]). Eco- 
logical communities comprising many species continue to 
exist under environmental perturbations due to complex 
"predator-prey" relationships maintaining their struc- 
ture. In particular, in microbial communities predator 
microorganisms act as homeostatic regulators to correct 
microbial imbalances and restor the balanced environ- 
ment |2j, density-dependent migration sustains the pop- 
ulations of birds, fishes, and insects [H, complex species 
interaction was experimentally revealed in the dynam- 
ics of ecosystems of several deserts in Arizona Q • There 
are a variety of social and economic systems such as large 
firms and enterprizes, product and service markets whose 
dynamics is governed by self-regulation. In these systems 
the behavioral standards, social norms, can be regarded 
also as the homeostasis characteristics (see, e.g., [EH, 0] ) . 
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When one of the homeostasis parameters deviates from 
its normal value and comes close to the destruction 
threshold the system has to respond to this event in or- 
der to prevent a further variation of this parameter. For 
a system with many elements its response to such per- 
turbations is usually implemented via reaction of a spe- 
cial subsystem governing the functioning of the system 
as a whole. In living tissue blood flow through the ves- 
sel network transporting various biochemical compounds 
governs the homeostasis of the cellular tissue. In a rela- 
tively large product market agents of various levels form 
a trade network and perturbations in the supply-demand 
equilibrium are damped, in principle, via the reaction of 
this trade network as a whole. We will call these subsys- 
tems life-support networks. 

A life-support network is to be of hierarchical struc- 
ture, because it not only controls the system homeostasis 
but also supplies the elements with the required "nutri- 
ents" and withdraws their "life activity" products. Typ- 
ically "nutrients" come into a system centrally from the 
external environment, so, before getting every element 
their flow branches many times. 

In the ideal case a life-support network has to operate 
in such a manner that every element be supplied with the 
amount of "nutrients" required for its current activity. 
However, for the product flow to very in a given branch 
of the terminal hierarchy level it has to change at all the 
levels. This flow redistribution in turn can induce the 
flow perturbations in the other terminal branches, neigh- 
boring and distant ones, thereby disturb the functioning 
of the corresponding elements and cause an erratic be- 
havior of the system as a whole (Fig. |T|) . 
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FIG. 1: (Color online) Illustration of the interference in the 
element functioning caused by the hierarchical organization 
of the life-support network; intensification of element 1 activ- 
ity requiring the "nutrient" flow to be increased disturbs the 
"nutrient" flow at elements 2. 



The existence of tremendous amount of natural sys- 
tems with hierarchical organization enables us to suppose 
that this destructive interference should be suppressed 
or, at least, depressed by the appropriate self-regulation 
processes. If there is no special unit governing the func- 
tioning of such system as a whole the self-regulation can 
be implemented only via a cooperative interaction be- 
tween all its elements. The latter feature allows us to 
call such regulation processes distributed self-regulation. 

Up to now the specific mechanisms by which the dis- 
tributed self-regulation arises are far from being under- 
stood well. In order to construct an appropriate model 
for the distributed self-regulation it is necessary to eluci- 
date two aspects. The first one is the self-processing of 
information. In a system without a governing center none 
of the elements forming the life-support network pos- 
sesses the complete information about the system state. 
Moreover, it really needs only a piece of information re- 
quired for the individual functioning. The second as- 
pect concerns the cooperative mechanism governing the 
product flow redistribution over the life-support network. 
None of these elements controls individually the product 
flow even through itself; responding to the appropriate 
information it can change only its own characteristics af- 
fecting the product flow. So the required redistribution 
of the product flow must be a cumulative effect of their 
individual contributions. 

Let us note some natural systems with the distributed 
self-regulation. First, it is living tissue, the main ob- 
ject of the present analysis, which will be discussed in 
detail in Sec. [Til Second, hierarchical organization is en- 
countered frequently in ecological communities based on 
trophic food pyramids Q (see also a review @). Their 
levels are made up of animals comparable in size and 
playing much the same role in the prey-predator rela- 
tionships. The linkages from small organisms generally 
vary over smaller scales. Larger animals that dominate 
these smaller organisms do so over larger scales. In 
such systems variations in species population disturbs lo- 



cally the prey-predator equilibrium and the induced self- 
regulation processes should sustain biodiversity [l(| EH ■ 
This concept is applicable also to socio-ecological systems 
comprising ecological communities, landscapes, and hu- 
man activities (see, e.g., fl2|). 

Third example is social networks describing many so- 
cial, ecological, and technological systems (for a review 
see [H, [lil JlaO . There are two categories of social 
networks [Ty] . One comprises the self-organized net- 
works, for instance, acquaintance webs or e-mail net- 
works, where hierarchy structures develop and evolve due 
to local agent interactions with the lack of global goals 
(see, e.g., [TtJ and references therein as well [THUD])- The 
other category is relevant to social systems with human 
collective efforts towards a given objective, e.g., large 
companies, distributed technology development, financial 
networks, and commodity markets or markets based on 
a certain raw material. The open-source communities 
16] as well as problem-solving networks aimed at prod- 
uct development [2(| also belong to the second category. 
The interplay between bottom-up decision making pe- 
riphery aimed at the overall goals of the community and 
top-down driven hierarchical organization can be imple- 
mented via response of such groups changing, for exam- 
ple, the number of programmers involved in their activity. 
A similar situation can be met in financial webs or mar- 
ket networks (see, e.g. @, Hl|). Hierarchical networks of 
administrative structures should have the capability of 
rearranging themselves to respond properly to failures, 
crises and disasters [Hj]. 

Papers [H, [H, [2f| [2(| [13] proposed several models for 
self-regulation of living tissue, a goods market based on 
a certain raw material, and ecological community. Their 
common feature is, first, the hierarchial self-processing 
of the information about the system state that is im- 
plemented via, for instance, conservation of blood flow 
and transported biochemical compounds or conservation 
of the raw material and money flow at the market net- 
work. Second, it is the individual response of, for ex- 
ample, blood vessels to biochemical compounds in blood 
flow going trough them or firms to disturbance in the lo- 
cal supply-demand equilibrium. It has been found that 
there exist special conditions under which self-regulation 
of such systems is perfect. Namely, if the functioning of 
an element has changed and it needs an increased amount 
of "nutrients" then the life-support network is able to do 
this in such a manner that the "nutrients" flow at the 
other elements be not disturbed at all. 

Unfortunately, the perfect self-regulation requires an 
ideal behavior of elements, so in real systems it cannot 
take place rigorously and should be regarded as the first 
approximation only. Since studying such systems at the 
initial stage it is worthwhile to single out some charac- 
teristic examples of these objects and to investigate them 
individually. The present paper considers living tissue in 
this way and analyze its response to local effects when 
the vessel properties are not ideal. 
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FIG. 2: Schematic illustration of the blood pressure distribu- 
tion over the vessels of systemic circulation and the roles that 
various groups of vessels play in blood flow distribution (after 
pfflj with some modification). 



FIG. 3: (Color online) Illustration of the vascular network 
layout made up of the counter-current pairs. The figure ex- 
hibits the rat renal vasculature reconstructed from computer 
tomography images (after [32], used with permission). 



II. PHYSIOLOGICAL BACKGROUND AND 
FUNDAMENTALS OF DISTRIBUTED 
SELF-REGULATION 

To justify the key points of the model to be developed 
let us note the relevant properties of the human physi- 
ology referring mainly textbook [Hj with respect to the 
general features. 



The architectonics of vessel network 

When an organ of human body is active and its tissue 
needs greatly increased supply of nutrients, blood flow 
can locally grow by 20 or 30 times. In contrast, the heart 
normally cannot increase its cardiac output more than 4 
to 7 times greater than when at rest. So, dilating or con- 
stricting, the blood vessels of each organ have to control 
locally blood flow at the required level. 

Figure [2] exhibits the distribution of blood pressure 
over the vessels of systemic circulation, which can be clas- 
sified into "conveying" and "delivering" types [30| . The 
large vessels, making up the "conveying" type, maintain 
a constant value of the mean blood pressure inside the 
large arteries and damp pressure pulsations via the spe- 
cial mechanism based on the baroreceptors located in the 
artery walls. As a result, at the entrance to the "deliv- 
ering" circulation system of individual organs the blood 
pressure can be considered as being fixed at the first ap- 
proximation if the whole organism is not under extreme 
conditions. 

Below the circulation system of an organ comprising 
small arteries and veins, arterioles and venules together 
with capillary network will be referred to as a microcir- 
culatory bed. As seen in Fig. [2] the main decrease in 
the blood pressure falls on the small arteries and arteri- 
oles and the blood pressure drop is rather uniformly dis- 



tributed over many vessels of different lengths. Thereby 
the net resistance of a microcirculatory bed to blood flow 
is determined by all its arteries and arterioles rather than 
vessels of one fixed length. 

In general, each artery entering an organ branches six 
to eight times before the arteries become small enough to 
be called arterioles, then the arterioles themselves branch 
two to five times before supplying blood to the capillar- 
ies. The arteries and arterioles are highly muscular and 
their lumen diameters can change manyfold. The branch- 
ing of these "delivering" vessels is rather symmetric; the 
daughter arteries of a mother artery are typically of the 
same size and form similar angles with the the parent 
vessel. 

Microcirculatory beds are typically arranged in a 
space-filling manner such that the terminal arterioles and 
venules be distributed in the tissue uniformly. Besides, 
rather often the layout of their arterial and venous net- 
works is mainly made up of the counter-current pairs 
consisting individually of an artery and vein located in 
space in close proximity to each other. Figure [3] illus- 
trates this counter-current vessel arrangement. Possible 
mechanisms of such pattern formation are discussed, in 
particular, in [3~il ]. 



The mechanisms of the tissue self-regulation 

To illustrate the state of art we discusses the tissue self- 
regulation mainly in the brain. In this organ the home- 
ostasis should be controlled by blood flow most effectively 
and its mechanisms have attracted much attention dur- 
ing the last decade (for a review see |33ll34l3al36l]). 

There is a lot of evidence that regulation of cerebral cir- 
culation is highly localized and cannot be governed solely 
by the products of tissue metabolism (see a review [37| 
and references therein). For example, during local acti- 
vation of a small region of the brain the corresponding 
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increase in the blood perfusion rate is caused by the di- 
lation of not only small arterioles in it but also relatively 
large pial arteries that supply directly the activated re- 
gion with blood [H, [H, Hcj ■ On one hand, exactly these 
arteries offer the greatest resistance to flow and, conse- 
quently, are the main site of flow control [4lj]. On the 
other hand, the dilated pial arteries lay outside the ac- 
tivated region and, thus, cannot be affected directly by 
biochemical products of cell activity in it. Thereby other 
mechanisms should also contribute to the blood flow reg- 
ulation. The main hypothesis about these mechanisms 
noted in many textbooks (e.g., [28|) is based on the re- 
lease of vasodilator substances such as nitric oxide (NO) 
by the cells of artery walls. Rapid blood flow through 
the arteries and arterioles causes shear stress on their 
cells, inducing the release of NO and finally the vessel 
dilation. The upstream propagation of this dilation is 
assumed to be conducted by the cell-to-cell interaction 
(through homoccllular gap junctions with the ions Ca 2+ 
playing the key role) or just the effect of the successive 
blood pressure redistribution (see, e.g., [33[). Yet these 
models are rather far from being verified well in vivo [36| , 
and other mechanisms are under discussion (see, in par- 
ticular, [H, |H, 0, HI, Hi|). Moreover, even the fact of 
observing in vivo upstream coupled waves of Ca 2+ and 
arteriolar dilation [47| poses a question as to how these 
waves propagating in both the directions along the arter- 
ies pass the points of vessel branching and transfer the 
information about the cellular tissue states. 

The upstream artery delation can be alternatively 
caused by artery-vein interaction; for a review see [481 ] 
and, in particular, paper [49| comparing sever al p ossible 
mechanisms with each other as well as work 50] study- 
ing this phenomenon in the brain. The main idea is that 
the downstream transport of biochemical substance with 
blood flow going through the vein network gives rise to 
an affective upstream processing of information about 
the cellular tissue state. The exchange of biochemical 
compounds between arteries and veins forming counter- 
current vessel pairs was assumed to arise via diffusion 
through the surrounding cellular tissue, being justified 
theoretically [HI [13]. 

Neural stimuli should also contribute to the local regu- 
lation of blood flow. In the brain the neural mechanisms 
of the blood flow regulation are due to the neurovascular 
units jH, [H, [H| . Roughly speaking, the neurovascular 
unit consists of vessels of microcirculation, interneurons, 
and astrocytes. Interneurons give rise to short pathways 
between between the tissue cells in the activated region 
and vessels supplying it with blood. Astrocytes together 
with interneurons affect the muscles of vessel wall causing 
dilation or contraction (for details see [5(| and references 
therein). 

The model to be developed actually combines the 
counter-current artery-vein interaction and the neural 
mechanism of local blood flow regulation (Fig. It 
assumes that the products of cell activity, for example, 
CO2, H + , and bioactive substances, including, the ions 
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FIG. 4: (Color online) Schematic illustration of physiological 
processes used in contracting the developed model. 



Ca 2+ , NO, etc. are accumulated in the cellular tissue 
and then with blood enter the vein bed. Below these 
biochemical substances are called activators. Receptors 
embedded into the vein walls detect the concentration O 
of the activators in blood flow going through the corre- 
sponding veins. Then the "reading" of the receptors are 
transferred by the neural mechanism to the arteries cou- 
pled with these veins, which determines the degree of the 
artery dilation or contraction. It should be noted that 
the substance diffusion between the nearest artery and 
vein can also imitate the action of these receptors. 

The next subsection demonstrates that the proposed 
model does give rise to the required upstream response 
of the arteries to disturbance in the cellular tissue home- 
ostasis. Moreover, if this artery response could be caused 
by the upstream propagation of some biochemical sub- 
stance along the arteries that obeys conservation at the 
branching points of the artery tree then the developed 
model could also allow for this mechanism after minor 
modification. 



The fundamentals of the distributed self-regulation 

Figure [5] illustrates the implementation of the self- 
processing of information about the states of living tis- 
sue. Let the activators be uniformly distributed in the 
domains Q\ and Q2 individually with concentrations Oi 
and O2, respectively. So the values Oi and O2 char- 
acterize the states of the living tissue in these domains 
and should determine the response of arteries a\ and 
ci2 directly supplying the living tissue domains Q\ and 
Q2 with blood, oxygen, and nutrients. Draining the 
given domains via the capillary network being in local 
quasiequilibrium with the surrounding tissue gives rise to 
the same values of the activator concentrations in veins 
V\ and V2 forming the counter-current pairs with the ar- 
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FIG. 5: (Color online) Illustration of the information self- 
processing caused by mass conservation in the blood flow 
through the venous bed. 



teries a\ and a^. The mass conservation at the shown 
branching points reads 



Jl2 
•^12©12 



(1) 

(2) 



where J%, J2, and J12 are the blood flow rates in these 
vessels. Whence it follows that the concentration of ac- 
tivators in the parent vein v\2 equal to 
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Ji6i + j 2 e 2 
J1 + J2 



(3) 



is the mean values of the activator concentration in the 
joint domain Q12 = Q\ 1J Q2 supplied with blood directly 
by the pair of the artery a±2 and the vein v\2- In this av- 
eraging the individual blood flow rates J\ , J2 play role of 
the weight coefficients. Thereby, the activator concentra- 
tion 0i2 is exactly the quantity characterizing the state 
of the living tissue domain Q12 as whole and required for 
the parent artery 012 to respond properly to variations 
in the state of living tissue in Q12. 

Generalizing this analysis and assuming the mass con- 
servation to hold at all the branching points of the vein 
bed we can write 



Jv — ^ ^ Jv c 

v c £v 



(4) 



where 0„ and J v are the concentration of activators and 
the blood flow rate at a given vein v and the sum runs 
over all the venules connected directly with the capillary 
network and originating from the vein v. For the sake 
of simplicity, the terms "arteriole" and "venule" will be 
used to refer to such smallest arteries and veins related 
connected directly with the capillary network, whereas 
the terms "artery" and "vein" will be corresponded to 
the other vessels of arterial and venous beds. Let Q„ c 
be the tissue region supplied with blood via the coupled 
arteriole and venule as a whole. All such domains will 
be called the elementary domains of living tissue. Then 
expression ^ gives us the mean value of the activator 
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FIG. 6: The resistance of an artery to blood flow through it 
vs the concentration of activators in the coupled vein. The 
model of perfect self-regulation assumes this dependence to be 
of the ideal form and the resulting dependence of the blood 
perfusion rate on the local value of activator concentration in 
the cellular tissue is shown in inset. 



concentration in the domain Q v = J2 V ev 2«c composed 
of all the elementary domains whose venules originate 
from the vein v. 

The receptors located in the vein walls detect the pres- 
ence of activators in the blood flow going through the cor- 
responding veins by measuring their concentration O v . 
Via the assumed communication between the artery and 
vein of a counter-current pair of vessels the receptor sig- 
nal is transmitted to the artery, causing its dilation or 
contraction to the degree determined by the value Q v . 
The hither the concentration, the stronger the artery di- 
lation should be. The effect of artery reaction will be 
measured in the dependence of its hydrodynamic resis- 
tance R on the activator concentration in the coupled 
vein. 

When one of the homeostasis parameters comes close 
to the boundary of the tolerance zone the tissue cells 
should generate an increased amount of the activators to 
signalize the system about the critical state, causing the 
blood flow rate to grow. For example, when the tissue 
temperature deviates from the normal value 36.6° C and 
gets the boundary about 44-45° C the blood flow rate has 
to increase essentially to prevent the further temperature 
growth. In the superficial tissue such temperature growth 
can force the blood flow rate to increase tenfold or higher 
[5^ . By the symbol A let us denote the concentration 
of activators corresponding to such critical situations. In 
other words, A specifies the threshold in the individual 
artery response; when the activator concentration in the 
coupled vein attains the value A, the artery has to dilate 
maximally. The general form of the adopted dependence 
of the artery resistance R to blood flow through it on the 
activator concentration in the coupled vein is depicted 
in Fig. [5] by solid line. Below as well as in Fig. [5] the 
activator concentration is measured in the dimensionless 
units, 8 = 0/A. 
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When the i?(#)-dependence is linear and can get zero 
value at 9 = 1 the vascular network is characterized by 
the perfect self-regulation [H, HH • Let us define the 
blood perfusion rate 77 at a point r of the living tissue as 
the ratio 
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venous bed 



rj(r) 



for Q Vc 3 r , 



(5) 




where V Vc is the volume of the elementary domain Q Vc 
containing the point r. In these terms the perfectness 
of the vascular network response is reduced to the state- 
ment that the blood perfusion rate r\ (r) depends only on 
the activator concentration 9(r) at the same point of cel- 
lular tissue. It is the result of cooperative interaction of 
all the arteries in governing the blood flow redistribution 
over the vascular network, which is illustrated in Fig. [7] 
To explain this fact let us assume that cellular tissue in 
the elementary domain Qo comes close to the tolerance 
zone and, as a result, the concentration of activators in it 
increases essentially. Then the activator concentration in 
all the veins forming up the path P" also grows, causing 
the dilation of all the arteries belong to the path P a . To 
analyze the effect of this dilation beyond the domain Qo 
we consider the induced variation of the blood perfusion 
rate at a point r ^ Qo- The widening of all the arteries 
located on the path P a before the branching point b° (i.e. 
between the entrance to the root artery and this branch- 
ing point) contributes individually to the increase of the 
blood perfusion rate at the point r because their dilation 
increases the total blood flow through the given vascular 
network. This effect is denoted by the sign "+" in Fig. [7] 
In contrast, the widening of the arteries located after the 
branching point b° have to decrease partially the blood 
perfusion rate r](r) because their dilation gives rise to a 
certain dominance of these arteries in the blood flow re- 
distribution after passing the branching point This 
effect is d esig nated by the sign "-" in Fig. [7] As proved 
in works [24l [25| all these contributions of the arteries 
compensate each other keeping up the value of the blood 
perfusion rate unchanged at the distant points. Natu- 
rally deviation of the _R(#)-dependence from the ideal 
form breaks this compensation. In this case the blood 
perfusion rate will depend on the activator concentration 
also at distant points of cellular tissue. 

Finalizing this subsection we also note the fact that in 
the case of perfect self-regulation the threshold A of the 
individual vessel response and a similar value Oth charac- 
terizing the response of the vascular network as a whole 
to increase in the activator concentration in cellular tissue 
coincide with each other. Indeed, the quantity 0th can 
be regarded as such a value of the activator concentration 
in cellular tissue that gives rise to a substantial increase 
in the blood flow resistance Z of the vascular network as 
a whole and, as result, to significant growth of the perfu- 
sion rate. In the case of the perfect self-regulation due to 
the feasibility for the vessel resistance to get zero value 



FIG. 7: (Color online) Schematic illustration of the vessel in- 
teraction giving rise the perfect self-regulation of living tissue 
when the i?(#)-dependence is of the ideal form. 

the blood perfusion rate rj depends on 9 as 



which allows us to set 9 t h = A. In the real situation such 
a singularity cannot arise and the inequality t h > A 
should hold. 



III. MODEL OF LIVING TISSUE 

The model for the living tissue response to local pertur- 
bations consists of three units for the vascular network 
architectonics, the regularities governing blood flow re- 
distribution, and the description of individual vessel be- 
havior. 



A. Vascular network architectonics 

The model for the vascular network architectonics as- 
sumes the following. 

First, the arterial and venous beds are individually of 
the tree form [58[. Second, the branching of vessels is 
symmetrical, i.e. at the branching points of the arterial 
bed the radii and lengths of daughter arteries are the 
same as well as are the angles formed by the daughter 
arteries with their mother artery. A similar statement 
holds also for the points of the vein merging. Third, the 
vessels of each level are distributed uniformly in the living 
tissue domain. This domain is considered to be a square 
or cub [59l ]. 

The first and second assumptions enable us to intro- 
duce the vessel hierarchy and to order all the arteries 
and veins according to their position in the hierarchi- 
cal structure. The arteries and veins reachable from the 
root artery or vein by passing the same number of the 
branching points (without return) belong to one level of 
the hierarchy. All the arteries or veins of one level are 
identical in properties. In what follows the root vessel 
will be labeled with index n = 0, the vessels of the next 
level are labeled with index n = 1, and so on. The third 
assumptions allows us to specify the vessel arrangement 
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FIG. 8: Dichotomic model of the vessel tree and its charac- 
teristic fragments of the embedding into 2D and 3D domains. 



in the living tissue domain using a self-similar embedding 
of the vessels. For example, the length of daughter arter- 
ies l n +\ is related to the length l n of the mother artery by 
the ratio l n /ln+i = / taking the same value / for all the 
levels. Rigorously speaking, this relationship holds only 
for the transitions between two or three successive levels 
for 2D and 3D models. However, in order to not overload 
description we will use this statement where it does not 
lead to misunderstanding. In particular, for the 2D and 
3D cases we set /213 = a/2 and Jzd = \ 2, respectively. 
Figure [H] depicts the characteristic fragment of the vascu- 
lar network that illustrates the embedding of vessels into 
the living tissue domain adopted in the present model. 
The total number of the hierarchy levels (not counting 
the root vessels) is divisible by 2 or 3 in 2D or 3D cases. 

Besides, for the sake of simplicity, the arterial and ve- 
nous beds are assumed to be the mirror images of each 
other with respect to all the properties, including the 
vessel response. This artificial assumption just simplifies 
the description of the vascular network response to per- 
turbations in the tissue homeostasis. Indeed, the vascular 
resistance of the real venous beds is small in comparison 
with that of the arterial beds. Therefore for a vascular 
network made up of the counter-current artery- vein pairs 
the blood flow patterns of the arterial and venous beds 
should be the mirror images of each other at the first 
approximation. So the latter assumption simply mimics 
this situation and enables us to speak only about veins 
and their dilation or contraction in the response to varia- 
tions in the activator concentration in blood flow through 
a given vein. Naturally, in order to get the right values 
of the blood flow rates in this model we have to increase 
the blood pressure drop across the vascular network twice 
or, what is the same, to fix formally the pressure inside 
the cellular tissue equal to the initial pressure at the en- 
trance to the root artery. Figure[9]illustrates the adopted 
model for the venous bed as well as the mechanism of its 
response to the variations in the activator concentration 
in the cellular tissue. 
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FIG. 9: (Color online) Vascular network model reduced to the 
venous bed and its reaction to the distribution of activators 
in the cellular tissue. Here, in particular, Q — Q Vc designates 
the elementary domain related to the venule v = v c - 

B. Governing equations for blood flow distribution 

Blood flow distribution over the vein network {v} 
(Fig. [5]) is described as follows. First, to blood flow with 
transported activators through every vein v we ascribe 
the blood flow rate J v and the activator concentration 
9 V . The state of this vein is quantified in terms of the 
vessel resistance R n (9 v ) depending on the activator con- 
centration 6 C in it, which reflects the vascular network 
response. Here the function R n (8) is assumed to the 
identical for all the vein of one level n. Second, to every 
branching point k we ascribe a certain blood pressure P. 
Then the blood flow distribution is governed by the sys- 
tem of equations including Ohm's law written for every 
vein v (Fig. [TO)) 

-Pin — -Pout = JvRn(9v) , (6) 

and conservation of blood and activator substance at ev- 
ery branching point k 

Jin,l + Jin,2 = Jout , (7) 
</in,l^in,l + </in,2#in,2 = */out#out ■ (8) 

This system of equations should be completed by the 
"boundary" conditions specifying the blood pressure at 
the entrance to the venules {v c } and the exit from the 
root vein 

Pin\v c =P, ^out|„=0-0 (9) 

and the equality of the activator concentration in every 
venule v c to the concentration of activators, or more rig- 
orously, its mean value in the corresponding elementary 
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FIG. 10: The main units of the mathematical description of 
blood flow distribution over the vein bed. 
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(10) 



The latter equality takes into account that blood flow 
through capillaries is rather slow and the capillary walls 
have numerous minute capillary pores permeable to wa- 
ter and other small molecular substances. This gives rise 
to local equilibrium between blood and the surrounding 
cellular tissue with respect to redistribution of biochem- 
ical compounds and makes the activator concentration a 
practically uniform held on scales about the size of ele- 
mentary domain. 

In the problem under consideration the distribution 
of activators is treated as a given beforehand held 9(r). 
Solving the system of equations I©-© subject to con- 
ditions ([9]) and (|10[) we can hnd the blood flow rate in 
the venules {v c }. The latter together with expression ([5|) 
gives us the desired dependence of the blood perfusion 
rate rj{9} on the held 9{r). 



C. Regularities of the vessel response 

To complete the developed model we should specify 
the dependence of the hydrodynamic resistance R V (9 V ) 
of every vein v on the activator concentration in it. 
According to the adopted assumptions the value 9 C = 1 
(O c = A in dimensional units) characterizes the criti- 
cal conditions when one of the homeostasis parameters 
comes close to the boundary of tolerance zone. To sur- 
vive the living tissue has to increase essentially the blood 
perfusion rate which is implemented via vessel dilation. 
Since all the arteries of a microcirculatory bed contribute 
equally to its vascular resistance each artery should ex- 
haust its ability to widen when the activator concentra- 
tion exceeds the critical value 9 C = 1. This allows us to 
adopt the following ansatz 



Rv{0v) = pn<f>(&v) ■ 



(11) 



where p n is the vessel resistance at 9 V — which de- 
pends only the number of the hierarchy level and 4>(6) 
is a certain function universal for all the vessel under 



FIG. 11: The expected dependence between the vessel dilation 
measured in units of <f> := R/Re=o and the receptor signal u(0) 
caused by activators of concentration 6. 



consideration. Naturally, by definition, at 9 — the 
equality 0(0) = 1 holds and for 9 3> 1 the inequality 
(f)(9) — > </>ii m <C 1 should be the case. 

In order to construct the function (f)(9) let us apply to 
the following speculations. There should be some thresh- 
old e -C 1 in the receptor ability to detect activators, 
meaning that activators cannot be detected if their con- 
centration 9 < e. So, to be specific, we suppose that the 
"signal" u(9) generated by the receptors is related to the 
activator concentration by the ansatz 



i{9) = V9 2 +e 2 -e. 



(12) 



Then the expected relationship between the function i 
quantifying the vessel response and the receptor "signal" 
u(9) should be of the form shown in Fig.[TlJ In particular, 
as u — > oo the value — > 0ii m and at u = u(9 = 1) the 
vessel dilation should get such a degree that the value 
be about its limit; we set 4>e=i = ^<piim- The following 
ansatz 



1 



20 



2 

lim 



"0=1 



l 



4> - <t>\im 



(13) 



has been chosen to specify this dependence. The combi- 
nation of expressions (fT2f and (fT3f yields the final ansatz 
for the vessel response function 



<km + U{ff) , [fam ~ U{0)] 



2*L, (14) 



where 
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mm 



20L 



9 2 



1 - 



'lim 



(15) 



Figure [12] visualizes the constructed function 4>(9) for 
various values of its parameters, including e = and 
<him = matching the ideal vessel response. In some 
sense, ansatz (|14l) can be justified applying directly to 
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FIG. 12: The vessel response function 4>(8) specified by 
ansatz (| 14[) for various values of its parameters. The upper 
fragment exhibits this function for <^ii m = 0.1 and e = 0.5, 
1.0, 10.0., the lower one does it for (j>y lm = 0.1, 0.2, 0.3 when 
the fixed value of e = 0.5. 



Fig. Q2] because the purpose of its construction is only 
to have the ability to analyze different cases of the living 
tissue response varying the parameters of an appropriate 
ansatz. 

The dependence of the vessel resistance on the hier- 
archy level, i.e. the value of p n vs the level number n 
is constructed using the statement about the equality of 
contributions to the vascular resistance of arterial beds 
(Sec. |nj from all the hierarchy levels. In other words, 
when the activator concentration is the same at all the 
points of the cellular tissue, in particular, 0(r) = 0, and, 
thus, all the vessels of one level have the same hydro- 
dynamic resistance, the blood pressure drop should be 
uniformly distributed over the vascular network. This 
enables us to write the /9„-dependence in the form 



Pn = (2C)>0 ■ 



(16) 



where the value po is related with the root vein and the 
parameter ( « 1 allows us to consider the cases where 
various groups of vessels dominate in the blood pres- 
sure distribution over the venous bed. In particular, for 
C < 1 the large veins contribute mainly to the vascular 
resistance of the vein bed, whereas, in the opposite case, 






FIG. 13: The structure of a vein branch B ViTl having the 
vein Vn as a root vessel (upper line) and the interconnection 
between the branches of neighboring levels (lower line). 



C > 1, the small veins are dominating. 



D. Numerical algorithm 

The developed model was analyzed numerically. So be- 
fore presenting the results the numerical algorithm used 
in solving equations dH])-© is briefly described. In addi- 
tion to veins we consider the collection of vein branches 
{B v }, each of which originates from one v of the veins 
and contains all its generations (Fig. The hierarchy 
of branches stems directly from the hierarchy of veins by 
ordering the branches according to the position of their 
root vessels in the vein tree. To every branch B v we 
ascribe the resistance to blood flow through it 



Zy 



Pn 



•J/J 



(17) 



just dividing the blood pressure drop across the branch, 
i.e. Pin,v (see Fig. [10]) by the blood flow rate J v at its 
root vein v n . To every branch B v we also ascribe the 
activator concentration 6 V in blood flow going through 
its root vein v. 

For the venules (veins of the last level) their branches 
{B Vc } consist only of these vessels. So their resistances 
)} can be found directly by virtue of con- 
dition (Tl0|) . To pass to the preceding level N — 1 we take 
into account the interconnection between the branches of 
the neighboring levels shown in Fig. [T3l Writing Ohm's 
law for these branches (see Exp. (fl7| ) and taking into 
account mass conservation at the branching point (see 
Exps. ([7]), ©) we, at the first step, can find the activa- 
tor concentration 6\2,n m the vein v\2, n from the formula 



'12, n — 



'1,71+1 



'2,n+l 



Tl+1 



J 2,n+1 



1 



1 



n+1 



J 2,n+1 



(18) 



if the values {9 n+ i} of activator concentration and the 
branch resistances {Z n+ i} are know at level (n+1), in 
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particular, we know the quantities 0i in +i, 02,n+i, 
and Z 2 .n+i ascribed to the branches B vi n+ i and B V2 n+1 . 
Then, at the next step, we calculate the resistance Z 12jl 
of the branch B Vl2 ra from the equality 
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Z12 
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n+1 



2, n+1 



(19) 



In this way all the characteristics of level (N — 1) can 
be found using these data of level N and so on up to the 
root vein of level 0. After reaching the root vein the total 
blood flow rate Jo = P/Zq is found. Then going in the 
opposite direction we calculate the blood flow rates { Ji} 
at level 1 using the equality 



Jl(2),n+1 — J 



1 



12, n 



Zi 



(2),n+l 



1 



1 



Z 



2,n+l 



(20) 



describing the blood flow split at branching points and 
so on until getting the venules. At the final step ex- 
pression (|5|) gives us the desired dependence of the blood 
perfusion rate on the distribution of activators in cellular 
tissue. 



IV. RESULTS OF SIMULATION 

The blood flow redistribution governed by equa- 
tions flS])-© subject to conditions © and (fTT)|) was stud- 
ied numerically using the algorithm described in the pre- 
vious section. Variations in the blood flow rate at every 
vein v as well as in the blood perfusion rate are analyzed 
in units of the corresponding quantities, /„ or 770, match- 
ing the absence of the activators, 9(r) = 0. These units 
allow us to reduce the number of the system parameters; 
only the total number of levels N, the parameter ( mea- 
suring the relative contribution from the lager and small 
vessels to the vascular network resistance, and the pa- 
rameters e and </>i; m of the individual vessel response are 
essential. 

The given section presents the results for the 2D sys- 
tem when the activator concentration 9(r) differs from 
zero inside m first elementary domains located at one of 
the square corner, namely, 



(r) 



0, r£ 



riQe 



Ui< 



Qi 



(21) 



Figure Q3] illustrates the analyzed case and shows the 
introduced ordering of the venules with index i. The 
elementary domains are also labeled with the index i of 
the corresponding venule. 

Figure [15] presents the obtained results of numerical 
simulation. The first column exhibits the dependence of 
the blood perfusion rate r]{9) on the activator concentra- 
tion 9 at the same point (through the same elementary 
domain) of living tissue. In obtaining these results acti- 
vators were assumed to be located only in the first ele- 
mentary domain Q\ (m — 1, except for the case shown 



in the left bottom frame) and the number of the hierar- 
chy levels was set equal to N = 6 (but the case shown in 
the left middle frame). It meets tenfold decrease in the 
vessels length when passing from the root vessel to the 
smallest arterioles or venules. 

The left upper frame depicts three curves correspond- 
ing to different values of the parameter £. When the ac- 
tivator concentration in the domain Q\ gets the thresh- 
old 9 C = 1 of the individual vessel response the venule 
joined with this domain exhausts its capacity to widen. 
However it is not enough for the blood flow rate to grow 
substantially. Indeed, the contribution of this venule to 
the vascular resistance is small even along the path Ve 
leading from the root vessel to the domain Qi on the 
vein tree (Fig. I14[) . So many veins along this path have 
to dilate in order for the blood perfusion rate to grow 
essentially in Q\. However, the larger a vein, the larger 
the tissue region drained through this vein as a whole 
and its size increases exponentially with the number of 
branching points separating the venule V\ and the given 
vein along Ve- As a result when the activator concen- 
tration in the domain Q\ is about 9 C = 1 the activator 
concentration in relatively large veins on the path Ve 
turns out to be rather small and these vessels cannot 
widen remarkably. Thereby the blood flow rate cannot 
also exhibit substantial increase (Fig. [T5|) . As the activa- 
tor concentration grows further the induced increase in 
the activator concentration along the path Ve gets the 
critical value 9 C = 1 in larger and larger veins, causing 
the growth of the blood perfusion rate in the domain Q± . 
When the main part of the veins belonging to Ve widen 
to the maximum the blood perfusion rate comes to the 
upper limit 77 max . 

The maximum ?7 max of the blood perfusion rate de- 
pends on the system parameter, in particular, on £. As it 
will be demonstrated below practically all the additional 
amount of blood entering the system due to the vessel 
dilation is directed to the domain Q\, So the higher the 
contribution of large veins, the higher the upper limit 
of the blood perfusion rate. These speculations are jus- 
tified well in Fig. [T5l the left upper frame. Because of 
this blood flow focusing along the path Ve the upper 
limit )j max of the blood perfusion rate can exceed sub- 
stantially the value l/0ii m which would be attained if 
the total blood pressure drop had fallen just on the last 
vessel supplying the region Qi with blood. 

A similar effect should be expected for the dependence 
of the blood perfusion rate on the number of hierarchy 
levels. The higher this number, the larger the upper limit 
?7max of the blood perfusion rate, which is demonstrated 
directly in Fig. [15l the left middle frame. When the re- 
gion of living tissue where activators are located becomes 
larger it effectively renormalizes the size of elementary 
domains and decreases the number of hierarchy levels. 
So the value of ?7 m ax has to drop with m increasing. The 
latter behavior is demonstrated in Fig. [TS] (left bottom 
frame), in particular, the value rj max goes down to (j)u m 
when activators spread over the whole region of the mi- 
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concentration gets the value 9 = 2 and 
the induced growth of blood perfusion 
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FIG. 14: (Color online) Example of the systems analyzed numerically. The region with a nonzero value of the activator 
concentration 6 is shadowed and the ordering of the vanules is illustrated for the system with the number of the hierarchy levels 
equal to N = 6. The right fragment visualizes the blood flow distribution in the domain of microcirculation bed based on the 
data obtained numerically and presented in Fig. 1151 



crocirculatory bed. 

The obtained results demonstrate us also the fact that 
when the vessel behavior is not ideal the threshold 9th in 
the vascular network response as a whole does not coin- 
cide with the critical value 9 C = 1 of the individual vessel 
reaction and can exceed it remarkably. Let us define the 
threshold 9th as the value at which the blood perfusion 
rate attains one half of its upper limit. Than this fact 
is clearly visible in Fig. [15] (left upper frame) and, for 
example, the estimate 9th ~ 2 holds for the case shown 
by curve 1. Naturally, the threshold 9 t h depends also on 
the parameters of the vascular network but not only on 
the characteristics of individual vessel behavior. 

The shown curves enable us to assume that due to 
the fractal structure of the vascular network there is a 
certain universal function r](9\(j), e) depending only on the 
parameters of individual vessel behavior, which describes 
the increase in the blood perfusion rate with the growth 
of the activator concentration until the latter becomes 
greater then the threshold 9 t h- The other characteristics 
of the vascular network, roughly speaking, affect only the 
upper limit ry max of blood perfusion rate. 

The right column depicts the distribution of the blood 
flow rates over the venules. Comparing the data shown 
hear and that of the left column we can declare that 
the main increase in the blood perfusion rate is located 
in the region of living tissue containing activators. Ex- 
actly these data enables us to assumed that the addi- 
tional amount of blood flow has to be directed mainly 



along the path Ve- However, the question about the in- 
terference of blood streams induced by distant regions of 
living tissue requires an individual investigation. Nev- 
ertheless, when the living tissue is excited only in one 
domain its response can treated as quasilocal, at least 
at the first approximation. We have used the quasilocal 
term because this response depends also on the size of 
the excited domain. 

As the nonlocal component of the ?7{#}-functional is 
concerned, the right column of Fig. [T5] demonstrates us 
that depending on the relation between the activator con- 
centration 9 in the elementary domain Q\ and the vas- 
cular network threshold 9th it changes its form. When 
9 < 9th the increase in the blood flow rate through Qi 
is caused by the blood flow decrease through the other 
elementary domains. In this case the required growth 
of the blood perfusion rate is due to the redistributed 
of blood flow over the vascular network without addi- 
tional amount of blood entering the microcirculatory bed. 
Since the number of the surrounding domains is rather 
large the decrease in the blood perfusion rate outside the 
excited region is not essential. When 9 > 9th the in- 
crease in the blood perfusion rate is mainly caused by 
additional portion of blood entering the microcirculatory 
bed. Therefore the blood perfusion rate is increased in 
the surrounding domains also. 

The resulting distribution of the blood perfusion rate 
in living tissue is shown in Fig. [T?] (right fragment) for 
the 2D model. As seen the distribution can be rather het- 
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FIG. 15: Results of numerical simulation. Blood perfusion rate r\ (in units of r?o) vs the concentration of activators in the same 
elementary domain Qi (left column) and the distribution of the blood flow rates (in units of Ji,N,e=o) over the venules, i.e. 
over the vessels of the last hierarchy level (right column). In obtaining these results the parameters e = 1.0 and 0i im = 0.1 were 
used, the used values of the other parameters are shown in figure. 



erogeneous because the elementary domains neighboring 
in the physical space can be supplied with blood through 
the veins (and arterioles) belonging to substantially dif- 
ferent branches. If activators are located in another ele- 
mentary domain the found distribution will also describe 
this case after a simple reordering of the vessels whereas 
in the physical space the blood perfusion pattern can 
become more heterogeneous then the shown one. The 
heterogeneity of the blood perfusion rate is caused by 
the presence of two topologies in the system, one is re- 
lated to physical space and the other corresponds to the 



vascular network architectonics. 

V. CONCLUSION 

The paper develops a model for the living tissue as an 
example of natural complex systems with active behav- 
ior. The necessity for such systems to keep their param- 
eters within tolerance zone makes their self-regulation 
crucial. Many of these systems like living tissue and 
consumer goods markets are organized hierarchically, i.e. 
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have highly branching networks supplying their elements 
with the required "nutrients". The absence of central 
governing units poses a question about the mechanism 
of self-regulation by which the "nutrients" are delivered 
to the points where it is necessary and the equilibrium 
at the other points is not disturbed. When an element 
requires additional amount of "nutrients" the supply net- 
work responding at all the levels should deliver them. 

Keeping in mind the previous papers devoted to this 
problem [23|, [24, 25, 26] we have developed a model for a 
microcirculation bed (a unit of regional circulation) em- 
bedded in 2D and 3D dimensional domains of cellular 
tissue in the self-filling manner. The cells of the tissue 
continuously require oxygen, various nutrient substances 
and the products of their live activity, for example, car- 
bon dioxide, should be withdrawn. When the intensity of 
cell activity grows the supply of these compounds has to 
be increased. Exactly this function of the vascular net- 
work is under investigation. We have assumed that the 
necessity for the increased blood supply is signalized by 
cells via the release of some activators, special biochem- 
ical compound. 

The artery and vein networks are considered to be of 
the tree form made up of the artery- vein counter-current 
pairs. Activators produced by cellular tissue are with- 
drawn by flood flow leaving the system via the venous 
bed. The considered mechanism of living tissue self- 
regulation is based on two effects. The first one is conser- 
vation of blood flow and the amount of activators at the 
branching points of the vein network. It gives rise to the 
information self-processing which is necessary for vessels 
to respond individually in an adequate way. The hydro- 
dynamic resistance of arteries is determined by the con- 
centration of activators in the blood flow going through 
the vein forming with a given artery a counter-current 
pair. For the sake of simplicity this model is reduced to 
the description of only the venous bed. The dependence 
of the vessel resistance on the activator concentration is 
constructed in a rather general way and deviates substan- 
tially from the ideal form analyzed previously [23l. [241. |25| . 

The following results have been obtained. First, it is 
shown that the threshold Oth hi the vascular network 
response to variations in the activator concentration as 



a whole can differ essentially from the critical value A 
of the individual vessel response and is determined not 
only by the parameters of the vessel behavior but also the 
characteristics of the vascular network architectonics. 

Second, due the fractal structure of the vascular net- 
work the dependence of the blood perfusion rate on the 
activator concentration in the cellular tissue is quasilo- 
cal, at least, at the first approximation. This means that 
the blood perfusion rate ?y(0) depends mainly by the ac- 
tivator concentration O taken at the same point of the 
living tissue. The prefix "quasi" takes into account that 
this dependence is affected also by the size of the excited 
living tissue region. The results of numerical simulation 
have demonstrated the fact that the ?7(8)-dependence is 
of a certain universal form determined only by the prop- 
erties of the individual vessel response for not too high 
values of the activator concentration, < A and only the 
upper limit (/ max of the blood perfusion rate is affected 
by the characteristics of the vascular network architec- 
tonics and the size of the excited region. It should be 
reminded that the blood perfusion rate under the nor- 
mal conditions is treated as a known value. This found 
quasilocality is explained by the fact that the main part 
of additional amount of blood going through the system 
is focused along a path on the vascular network connect- 
ing the excited region and the root vessels. 

Third, due to the presence of two topologies in living 
tissue, the point proximity in the physical space and the 
vessel proximity with respect to their position in the vas- 
cular network the nonlocal component of the functional 
?7{0}-dependence can be rather heterogeneous in space 
because neighboring points of living tissue can be sup- 
plied with blood by vessels belonging to vessel branches 
distant from each other in the vascular network architec- 
tonics. 
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